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Numerical simulations of the complex cubic-quintic Ginzburg-Landau equation (CCQGLE), a 
canonical equation governing the weakly nonlinear behavior of dissipative systems in a wide variety 
y—{ of disciplines, reveal five entirely novel classes of pulse or solitary waves solutions, viz. pulsating, 

creeping, snaking, erupting, and chaotical solitons [9]. Here, we develop a theoretical framework 
^SJ for analyzing the full spatio-temporal structure of one class of dissipative solution (snaking soliton) 

, ^ of the CCQGLE using the variational approximation technique and the dynamical systems theory. 

^ ^ The qualitative behavior of the snaking soliton is investigated using the numerical simulations of (a) 

the full nonlinear complex partial differential equation and (b) a system of three ordinary differential 
^ equations resulting from the variational approximation. 

oo 

^ I. INTRODUCTION 

The complex cubic-quintic Ginzburg-Landau equation (CCQGLE) is one of the most widely studied non- 
• 1-H linear equations. In fluid mechanics, it is also often referred to as the Newell- Whitehead equation after the 

(— I authors who derived it in the context of Benard convection. Many basic properties of the equation and its 

'— ' solutions are reviewed in [2 , together with applications to a vast variety of phenomena including nonlinear 

^ waves, second-order phase transitions, superconductivity, superfluidity, Bose-Einstein condensation, liquid 

crystals and string theory. 

Due to the highly nonlinear nature of the CCQGLE, analytic solutions have been limited to special cases. 
Thus, its numerical approximations have become valuable tools in order to better understand the theoretical 
T-H insight into the problem. Numerical simulations in identify five entirely novel classes of pulse or solitary 

waves solutions, viz. pulsating, creeping, snaking, erupting, and chaotical solitons of the CCQGLE. In 
contrast to the regular solitary waves investigated in numerous integrable and non-integrable systems over 
the last three decades, these dissipative solitons are not stationary in time. Rather, they are spatially 
confined pulse-type structures whose envelopes exhibit complicated temporal dynamics. 

In this paper we describe a theoretical framework for analyzing the full spatiotemporal structure of one 
^ class of solitary waves (snaking soliton) of the CCQGLE, and present numerical simulations resulting from 

the full nonlinear partial differential equation, as well as systems of DDEs from the variational approximation 
technique. 

First, we develop and discuss a variational formalism within which to explore the various classes of dissipa- 
tive solitons. We find a suitable Lagrangian based on ansatz. Then the resulting Euler-Lagrange equations 
are treated in a completely novel way. Rather than considering on the stable fixed points which correspond 
to the well-known stationary solitons or plain pulses, we use dynamical systems theory to focus on more 
complex attractors viz. periodic, quasiperiodic, and chaotic ones. Periodic evolution of the trial function 
parameters on a stable periodic attractor would yield solitons whose amplitudes are non-stationary or time 
dependent. In particular, pulsating, snaking (and less easily, creeping) dissipative solitons may be treated 
using stable periodic attractors of various trial function parameters. Chaotic evolution of the trial function 
parameters would yield to chaotic solitary waves. In the language of the Los Alamos school, the fully spa- 
tiotemporal approach followed here may be said to be the "collective coordinates" formulation. In other 
words, we consider a pulse or solitary wave at any time as a coherent collective entity (or coordinate). This 
solitary wave is then temporally modulated. The spatial approach proposed, and explored, in this paper 
is the variational method. However, the method is very significantly and non-trivially generalized from all 
earlier applications to deal with our novel classes of dissipative solitary waves. 

Next, we solve (a) the initial boundary value problem of the CCQGLE and (b) the initial value problem of 
the Euler-Lagrange system resulting from the variational approximation numerically. The numerical scheme 



X 



for the PDE is based on a Finite Difference discretization of the partial differential equation with a lo"^ 
storage varient explicit Runge-Kutta third order method. The Euler-Lagrange system of ODEs is solved 
numerically using the MATLAB's built-in function ode45. 

The remainder of this paper is organized as follows. Section §2 describes briefly the CCQGLE as well as 
the classes of its soliton solutions. An ODE system resulting from the generalized variational formulation of 
the CCQGLE and Hopf bifurcation theory are also presented here. In section §3 numerical methods for the 
IVP and BVP are discussed. Sections §4 and §5 deal with the simulation results from previous section using 
the PDE and the variational formulation. Finally, in section §5 we summarize the results and conclude the 
paper outlining some further research directions. 



II. DESCRIPTION OF CCQGLE 



An important element in the long time dynamics of pattern forming systems is a class of solutions which 
we call "coherent structures" or solitons. These structure could be a profile of light intensity, temperature, 
magnetic field, etc. A dissipative soliton is localized and exists for an extended period of time. Its parts 
are experiencing gain/loss of energy or mass with the medium. Thus, energy and matter can flow into the 
system through its boundaries. As long as the parameters in the system stay constant, the structure that 
could evolve by changing shape exists indefinitely in time. The structure disappears when the source is 
switched off, or if the parameters are moved outside the range of existence of the soliton solutions. Since 
these dissipative systems include energy exchange with external sources, they are no longer Hamiltonian, 
and the solitons solutions in these systems are also qualitatively different from those in Hamiltonian systems. 
In Hamiltonian systems, soliton solutions appear as a result of balance between diffraction (dispersion) and 
nonlinearity. Diffraction spreads a beam while nonlinearity will focus it and make it narrower. The balance 
between the two results in stationary solitary wave solutions, which usually form a one parameter family. In 
dissipative systems with gain and loss, in order to have stationary solutions, the gain and loss must be also 
balanced. This additional balance results in solutions which are fixed. Then the shape, amplitude and the 
width are all completely fixed by the parameters of the dissipative equation. However, the solitons, when 
they exist, can again be considered as "modes" of dissipative systems just as for non-dissipative ones. 

Recent perturbative treatments based on expressions about the nonlinear Schrodinger equations are gen- 
eralized to perturbations of the cubic-quintic and derivative Schrodinger equations. The cubic Ginzburg- 
Landau admits a selected range of exact soliton solutions T]. These exist when certain relations between 
parameters are satisfied. However, this certainly does not imply that the equations are integrable. Dissipa- 
tive systems, as it is the case with Burger's equation are integrable. In reality, general dissipative nonlinear 
PDEs cannot be reduced to linear equations in any known way, therefore an insight to the type of solutions 
that the equation has may be based on numerical simulations. 

The simplest mode that can account for this type of behavior is the Ginzburg-Landau model, with the 
corresponding cubic-quintic equations in the form [13] 

dtA = eA+ (fci + ici)dlA - (63 - jc3)|ApA - (65 - ic5)\A\''A. (II.l) 

The present study will confine itself to spatially infinite systems in one dimension and will focus primarily 
on a spatio-temporal behavior of dissipative solitons. One way to approach this problem and to arrive at 



(II.l I is to start with the well known modified Schrodinger equation 

dtA = ic^dlA + zcalApA + ic5\A\'A + a,[(so + S2\A\'')AI {11.2) 



and to perturb with a dissipative term in special parameter regimes, with terms on the rhs. of (II. 2) of 
the form f{A) ^tA + hd^A - b3\A\'^A - b5\A\^A. 

The unperturbed system ( |II.2[ ) leads to an integrable dynamical system whose orbits can be calculated 
analytically [J- 



The interpretation of the system's parameters in (II.l) depends on the particular field of work. For 



example, in optics, t is the propagation distance or the cavity number, x is the transverse variable, the 
angular spectral gain or loss is identified by 61, ci is the second-order diffraction or linear dispersion, 63 
is the nonlinear gain or 2-photon absorption if negative, C3 is the nonlinear dispersion, e is the difference 
between linear gain and loss, 65 represents the saturation of the nonlinear gain, and C5 the saturation of the 
nonlinear refractive index. In physical problems, the quintic nonlinearity is even of higher importance than 
the cubic one, as it is responsible for stability of localized solutions. 



A special case of (II. 1 1 is the nonlinear Schrodinger equation (e = &i = &3 = &5 = C5 = 0, ci = 1) 

dtA = id^A + ic3\A\^A (II.3) 

which is both Hamiltonian and integrable. Its extension, the quintic-cubic Scrodinger equation {e — bi = 

b3 = b5= 0, ci = 1) 

dtA = idlA + icM?A + ic^\AfA (II.4) 
is Hamiltonian but non-integrable. Other interesting cases are the derivative nonlinear Schrodinger equation 

dtA = idlA + s^d^A + s2d^(\A\^A) (II.5) 
and the combinations of last two, which is the quintic derivative Scrodinger equation 

dtA = idlA + s^d^A + s-2.d^{\A\^ A) + ic-i\Af A + ic^\A\^A. (II.6) 

A. Euler-Lagrange Equations 

Employing the generalized variational formulations and proceeding as in [8J, the Lagrangian for the CC- 



QGLE ( |II.1| may be written as [TO] 

L^r* {dtA -eA- {b, + ic^)dlA + (63 - ic3)\A\^ A + (65 - ic^)\A\^A\ 
+ r[5tA* - eA* - (&i - zcOa^A* + (63 + «c3)|A|2A* + (65 + ic^)\AfA*\ . (II.7) 

Here r is the usual auxiliary equation employed in [8] and it satisfies a perturbative evolution equation 
dual to the CCQGLE with all non-Hamiltonian terms reversed in sign. Choosing the single-humped trial 
functions of the form: 

A(x,t) = Ai(<)e-"i(*)'["-^i(*)l'e^"i(*)e*'^(*)'"' (II.8) 

where the A\(€) is the amplitude, the cri(t)'s are the inverse widths, (/)i(t)'s are the positions (with ^i(i) 
the speed), e*'''^*-' ^ represents the chirp, Q!i(t)'s are the phases of the solitons and are all allowed to vary 



arbitrarily in time. For now, the chirp terms are omitted for simplicity. Substituting (II.8|/(II.9) in (II. 7) 



the effective or averaged Lagrangian is 



L 



EFF = 



Ldx = 2\/Tr< — 



[cTl(i)2+a2(i)2]5 



eAi(i) cos[ai(i) - a2{t)] 



_ 3<7i(t)^,T2(') [■/'l(*)-*2(*) 
e 3,Ti(t)^+,T2(')^ 



3fTi(i)2+fT2(i)2 



+ 



5gi(t)"'g2(t)^[^l(t)-.»2(*)l- 
5^l(t)^ + ^2(t)^ 



5<Ti(i)2+cT2(i)2 

^l(t)^<^2(')''['»l(t)--»2(*)l^ 



63 cos[Q;i(t) - a2(t)] + C3 sin[ai(i) - Q;2(t)] 



65 cos[Q;i(t) - a2(t)] + C5 sin[ai(t) - a2{t)] 



cos[ai(i) - a2(t)][cri(i)^ + cr2{t)^]'^Ai{t) 



ai(t)2 + a2(t)2 

+ ^ - 2(Ti(i)V2(t)^ 61 cos[ai(^) - a2(t)] - ci sin[ai(i) - a2(t)] - cr2(t)^ 

+ ai(t)2[-l + 2a2{tf[Mt) - Mt)?]] - sin[ai(t) - a2{t)][aiit)^ + (XaW']' 

- ai(t)ai(i) cos[ai(t) - a2{t)] [ai{t)^ + <J2{t)^ + 2<j2{t)^[Mt) ' Mt)? 

- 24>i{t)ni{tfa2{tf[Mt) - (t)2{t)][ai(tf + n2{tf] cos[ai{t) - a2(t)] 



(11.10) 



Since (11.10) reveals that only the relative phase alt) — ai{t) — a2{t), and relative velocity = 4'i{t) — (l)2{t) 
of A{x, t) and r{x, t) are relevant, we henceforth rescale them according to 



ai(t) = a(t),a2W =0; (/.i(t) = 02 W = 0. 
Also, for algebraic tractability, we have found it necessary to assume 

02{t)=ai{t)=G{t). 



(11.11) 



(11.12) 



While this ties the widths of the A{x^ t) and r{x, t) fields together, the loss of generality is acceptable since 
the field r{x,t) has no real physical significance. Moreover, for the snaking soliton solutions, we choose 
a{t) — 2/(p{t). By using all these assumptions, the trial functions (II.8|-(II.9) become 



A{x,t) = Ai{t)e 



4 J.2 

r{x,t) = e ^c)^ . 



(11.13) 
(11.14) 



Substituting the last two equations into (II. 10 1, the effective Lagrangian becomes 

Leff = i^—T 6e3^i(t)30(t)2(63cosa(i) + C3sina(i)) 

12e 3 

cosa(t) + C5 sina(t)) 
- 3V2e-i ^ - 2Ai{t) sma{t){ - 12ci + (t>{tfa{t)) 

+ cosa{t){ - 2(j)'^{t)a(t) + Ai{t){24bi + 2e(j)'^{t) + 3(/.(<)(/>(i))) 



(11.15) 



We are left with three parameters Ai{t), (f){t) and a{t) in Leff- Varying these parameters by using calculus 
of variations, we obtain 

OLeff d /9L 



d*{t) dtV dMt) 

where ★ refers to Ai, a, or 0. Solving for *(t), we obtain the following system of three ordinary differential 
equations, 

ii(i) = /i[^i(i),a(t),(/.(t)] 
a{t)^ f2[Ai{t),a{t),(l){t)] 

Ht)^f3[Ai{t),a{t),^it)] (11.17) 
where the right hand side /i, /2, /s of ( II.17[ ) are given by 

/i = ^Ai(t)secQ;(i)(3456cie*/^sinQ;(t)(-5iCosa(t) + cisina(t)) 
5 

+Ai{t)^(j){t)'^{2V6e^^^Ai{t)^{5b3b5 + 17C3C5 + (86365 - 19C3C5) cos(2a(i)) 
+(-26365 + 2C3C5) cos(4a(i)) + (65C3 + 63C5)(11 sin(2a(t)) - 2 sin(4a(t))))(/>(t)^ 
+4Ai(t)6(46^ + 13cl + (26^ - Ibcl) cos(2a(t)) + 2(-6^ + cj) cos(4a(t)) 
+I765C5 sin(2a(i)) - 465C5 sin(4a(t)))(^(t)2 + Se'^/^Ai{tf{3bl + llc^ 
+2(6^ - 6cl) cos(2a(t)) + (-6^ + cj) cos{4a{t)) - 263C3(-7sin(2a(t)) 
+ sin(4a(t))))0(i)2 + 6\/2e5/3(66i63 - I8C1C3 + 6(6163 + ciC3)(2 cos(2a(i)) 
+ cos(4Q!(i))) — 48(6301 — 61C3) cosa(i)'^ sina(t) — 2ecos a{t){b3 cos a{t) 
+3C3 sina(t))(/)(i)2) + AV3e^^^Ai{t)^{12bib5 - SGciCg 

+12(6165 + ciC5)(2cos(2a(t)) + cos(4a(t))) - 96(6501 - 61C5) cos a(t)^ sin a(i) 
-2e cos a(t)(265 cos a{t) + 5c5 sin a{t))(j)(tf)))e-^/^ /(j)itf /{UAcie'^^^ sin a(t) 
+Ai(0'(-3y2ei/3(-363 cosQ!(t) + 63 cos(3Q(i)) — IIC3 sma{t) + C3 sin(3Q;(t))) 
+4^3^1(^)2(465 cos a{t) — 265 cos(3a(t)) + 13c5 sin a(t) 

-2c5sm{3a{t))))(j){tf) (11.18) 



/2 = — (-72e''/^(6i - citana(i)) + (-6e'*/^e + 9\/2ei/3Ai(<)2(63 + C3 tana(i)) 
15 

+10V3Ai{t)'^{b5 + C5tana{t)))(l){tf)e-'^^^/(l){t) (11.19) 



/3 = -(-144cie^/^ecosa(i) - V6Ai(t)^((65C3 - 363C5)cosa(0 
5 

+ (6503 + 63C5) cos(3q!(<)) - 2(6365 + C3C5 + (6365 - C3C5) cos(2a(t))) sin a(t))(/)(t)2 

+6y2e'*/^^i (t)2 (12(6301 +361C3) cos a(i) + 12(6301 - 61C3) cos(3a(i)) 

+12(6163 - 3ciC3)sinQ;(i) + 12(6163 + C1C3) sin(3Q;(i)) + e cos a(t)(c3(-2 

+ cos(2a(t))) - 63 sin(2a(t)))(/)(i)2) +8V3eAi(i)'*(366iC5 cos a{t) 

+18(6501 — 61 C5) cos(3Q!(i)) + 18(6165 — 3C1C5) sina(i) + 18(6165 + C1C5) sin(3Q:(t)) 

+ecosa(t)(-3c5 + 2c5 cos(2a(i)) - 265 sin{2a{t)))(l){tf))e~^/{UAcie'^/^ sina{t) 

+Ai{t)^{-AV2e^^^{-3b3 cosa{t) + 63 cos(3Q!(t)) — IIC3 sinQ!(t) 

+C3sin(3a(t))) +4%/3Ai(t)2(465 cosQ!(i) — 265 cos(3a(t)) + 13c5 sina(i) 

-2c5sm{3a{t))))(j){tf) (11.20) 



B. Hopf Bifurcations 
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The Euler-Lagrange equations (11.17) are treated in a completely novel way. Rather than considering 



the stable fixed points which correspond to the well-known stationary solitons or plain pulses, we use Hopf 
bifurcation theory to focus on periodic attractors. Periodic evolution of the trial function parameters on 
stable periodic attractors yields the pulsating soliton whose amplitude is non-stationary or time dependent. 

We derive the conditions for the temporal Hopf bifurcations of the fixed points. The conditions for super- 
critical temporal Hopf bifurcations, leading to stable periodic orbits of Ai{t), 4'{t), and a{t) can be evaluated 
using the method of Multiple Scales, as in [TT]. These are the conditions or parameter regimes which exhibit 
stable periodic oscillations, and hence stable pulsating solitons will exist within our variational formulation. 
Note that, it is easy to verify numerically, periodic oscillations of Ai{t), 4i(t), and a{t), correspond to a 
spatiotemporal pulsatin g solit on structure of the |A(a;,t)| given by ( |H.8[ ). 

The fixed points of (H.17) are given by a complicated system of transcendental equations. These are 



solved numerically to obtain results for each particular case. 

For a typical fixed point, the characteristic polynomial of the Jacobian matrix of a fixed point of (H.17) 
may be expressed as 



+ Si\^ + S2X + S3 = 



(11.21) 



where Si with i = 1...3 depend on the system parameters and the fixed points. Since these are extremely 
involved, we omit the actual expressions, and evaluate them numerically where needed. 

To be a stable fixed point within the linearized analysis, all the eigenvalues must have negative real parts. 
Using the Routh-Hurwitz criterion, the necessary and sufficient conditions for (11.21) to have -Re(Ai.2,3) < 
are: 



Si > 0, ^3 > 0, S1S2 -S3>0. 



(11.22) 



On the contrary, one may have the onset of instability of the plane wave solution occurring in one of the 
two ways. In the first, one root of (H.17) (or one eigenvalue of the Jacobian) becomes non-hyperbolic by 
going through zero for 



<5,3 - 0. 



(H.23) 



Equation (11.23) is thus the condition for the onset of "static" instability of the plane wave. Whether this 



bifurcation is a pitchfork or transcritical one, and its subcritical or supercritical nature, may be readily 



determined by deriving an appropriate canonical system in the vicinity of (H.23) using any of a variety of 
normal form or perturbation methods. 

One may also have the onset of dynamic instability ( "fiutter" in the language of Applied Mechanics) when 
a pair of eigenvalues of the Jacobian become purely imaginary. The consequent Hopf bifurcation at 



S1S2 







(H.24) 



leads to the onset of periodic solutions of (11.17) (dynamic instability or "flutter"). 



III. NUMERICAL METHODS 



A. The Initial Boundary Value Problem of the CCQGLE 



To obtain the snaking soliton solutions that are periodic and exponentially decaying at infinity, we develop 

A{b, •), and A{0, ■) — Ao{-) based on the finite 



numerical schemes for the IB VP of (II.l) with A{a, •) 
difference method as described below. 

Let h = {b — a)/n be the grid spacing and Xj — 
fine Aj(t) as an approximation to A{xj,t), j = 0,1,-' 
(Aoixi),--- ,Aoixn-i)V. Let F{Aj{t)) = eA,(t) - (63 



J- + jh, j 0, 1, • • • 
,n, A{t) = (Aiit),-- 
tC3)\A,{t)\^A,{t)-{b, 



be the grid points. De- 
,A„_i(i))^ and Ao(a;) = 
-ic5)\Aj{t)\^Aj{t). Using 



the central difference approximation to d^A{xj ,t), we write the semi-discretization of the IBVP of (II.l) as 



the following system of ODEs 
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A(t) =SA(t) +F(A(i)), Vi > 0, (III.l) 
A(0) = Ao, (III.2) 



where F{A{t)) = {F{Ai{t)), ■■■ , F{An{t)))'^ . The upper dot in (III.l) indicates derivative with respect 



to t and B is the finite difference matrix. To construct an integration scheme to solve the ODE system 
(III.l )-( III.2[ ), let tn+i = tn + At, and let A" denote the value of the variable A at time t„. Employing a 



low storage varient third-order Runge-Kutta scheme [Uj, we write the fully discrete system as 

Qi=AtG(A"), Ai-A" + iQi, 

Q2 = AiG(Ai)- |Qi, A2=Ai + ||Q2, (III.3) 
Q3 = AtG(A2) - |iQ2, A"+i = A2 + ^Q3, 

where G(A"') — SA" + F(A"). The numerical code is parallelized for distributed memory clusters of 
processors or heterogeneous networked computers using the MPI (Message Passing Interface) library and 
implemented in FORTRAN. 

Computations were performed on a Linux cluster (zeus.db.erau.edu: 256 Intel Xeon 3.2GIIz 1024 KB 
cache 4GB with Myrinet MX, GNU Linux) at Embry-Riddle Aeronautical University. 

B. The Initial Value Problem of Euler-Lagrange System 



The system of ODEs (11.17) resulting from the variational approximation to the PDE is solved numeri- 



cally using the MATLAB's built-in function ode45 (adaptive fourth- and fifth-order Ru nge-K utta-Fehlberg 



method). We use equilibrium solution Ai{0) — A^, a{0) — a°, (j){0) — 0° of the system (11.17) as the initial 
conditions. The numerical solutions A", a", giving the approximations to Ai{tn), a(i„), 4>{tn) are used 
to evaluate A{xj,tn), approximation to the snaking soliton solution to the CCQGLE, in the spatial grids xj 
at time tn using the following formula 

where a = 2 /</>". 

IV. RESULTS AND DISCUSSION 

A. Simulations of Snake Solitons using the PDE 

The numerical scheme for the PDE described in section §3 is implemented with the following initial 
amplitude profile 

A{x, 0) = Aie-'''(^-^)'e^". (IV.l) 
We use Ai — 0.583236, (p = 1.05969, a = 0.185515 and a = 1.8873 as the typical values of the parameters in 



(IV.l ). These are the values of the fixed point as we will see in the next section. For the system parameters 
of the PDE ( II.1[ ), initially we use the values listed in Table 1 (see the row for snaking soliton). In [9] five 



novel classes of dissipative soliton solutions viz. pulsating, creeping, snaking, exploding and chaotical were 



obtained by numerical simulations of the CCQGLE (II.l ). Here we present a theoretical formulation to one 



class of the solutions, the snake, we perform independent simulations on the full PDE from the ones shown 
previously in [Tj, and we also compare the result with the simulations from the variational approximation. 
First, for the simulations on the PDE, we fix a set of system parameters e, &i, b^, 65, Ci, C3, C5 for the 
snake soliton from the Table 1. Then, we study the qualitative behavior of the snaking soliton by varying 
one parameter at a time. Then, we present the spatio-temporal structure of the snaking soliton in detail 
using the numerical simulations of the variational approximation to CCQGLE in the next section. The 
spatiotemporal structure of the solitons obtained from the simulations of the partial differential equation 
(CCQGLE) is shown in Figs. 1-4. In each of the solitons presented in Figs. 1-4 one of the parameters is 




FIG. 1: Numerical simulation of snaking soliton using the PDE; unperturbed case 



TABLE L Parameters of the CCQGLE 



Solitons 


e 61 ci 63 C3 65 C5 


Pulsating 


-0.1 0.080 0.5 -0.660 1 0.10 -0.100 


Creeping 


-0.1 0.101 0.5 -1.300 1 0.30 -0.101 


Snaking 


-0.1 0.080 0.5 -0.835 1 0.11 -0.080 


Exploding 


-0.1 0.125 0.5 -1.000 1 0.10 -0.600 


Cliaotical 


-0.1 0.125 0.5 -0.300 1 0.10 -1.000 



perturbed while all the other remaining parameters are chosen from the Table 1 (see the row for the snake 
soliton, which is the unperturbed case). Thus, only the perturbed parameter is indicated in the captions to 
the figures while the other ones stay unchanged. 




FIG. 2: Perturbed snaking soliton using the PDE: Left 63 = -0.8, Right 65 = 0.15 



V. SNAKE SOLITONS USING VARIATIONAL APPROXIMATION 



In our previous work [TU], we looked at stationary soliton solutions of the CCQGLE in the parameter 
regimes where they exist. However, these are not the only possible type of solutions. Pulsating solutions 
are another example of localized structures. They arise naturally from the stationary ones when the latter 
become unstable. Pulsating solitons might have several frequencies in their motions and the solutions will 
be quasiperiodic. A relative simpler case is when the motion has two frequencies. Then, we will have a 
pulsating soliton |11| . which instead of having a zero velocity will move back and forward around a fixed 
point. Obviously, there are two frequencies involved in this motion which usually are incommensurate. We 
call this type of solitons with more than one frequency a snaking soliton. To capture this, we will use 
dynamical systems theory to construct solitons with quasiperiodic behavior. 

Stable fixed points of the Euler-Lagrange system (11.17 1 corresponds to the well-known stationary solitons 
or plain pulses. The periodic evolution of Ai{t), (l){t), a{t) on stable periodic attractors yields solitons whose 



FIG. 3: Perturbed snaking soliton using the PDE: Left Ci — 0.55, Right C3 — 1.135 





FIG. 4: Perturbed snaking soliton using the PDE: Left cs 



-0.06, Right e ■ 



-0.08 



amplitude is non-stationary. To derive the conditions for occurrence of stable periodic orbits of Ai{t), (j)(t), 
and a{t), we proceed as follows. 

First, we fix a set of system parameters 61 = 0.08, 65 = 0.11, ci = 0.5, C3 = 1, — —0.08 as in the 
case of the full PDE presented in the Table 1. Then, we solve numerically the system of transcendental 
equations (ILIT), which are the equations of the fixed points. By the Ruth-Hurwitz conditions, the Hopf 
curve is defined as S1S2 — =0. This condition, along with the equations of the fixed points leads to 
onset of periodic solutions of (11.17) as we will see next. On the Hopf bifurcation curve we obtain that 
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FIG. 5: Numerical simulation using the variational approximation: periodic time series of the amplitude Ai{t), the 
phase a{t) and the position (j){t) on the Hopf curve (Left), spatiotemporal structure of the amplitude |A(x,t)| = 
^„^-a^(^^-4,"y on the Hopf curve (Right). 



bsHopf = -1.89646, and eRopf = -0.297393, while the fixed points are Ai{0) = 0.583236, 0(0) = 1.05969, 
and a(0) = 0.185515. For these values of b^Hopf and SHopf, which are two free bifurcation parameters, 
we integrate numerically the systems of 3 ODEs (II. 17 1, using as initial conditions the three values of the 
fixed points, to find the trial function's periodic functions. Hopf bifurcations occur in this system leading 
to periodic orbits. These are the conditions which exhibit periodic oscillations, and hence stable snaking 
solitons will exist within our variational formulation. 



Next, we plot the time series for the ampHtude Ai{t), the phase a{t) and the position (p{t) on the Hopi 
curve and the results are presented in Fig. 5 (Left). As expected, we noticed that the amplitude is small, 
since the maximum height of the snaking soliton is proportional to the square root of the distance from the 
Hopf curve. The corresponding spatio-temporal structure of the ansatz |A(a;,i)| = A'^e~'^ i^j~4' ) on the 
Hopf curve is shown in Fig. 5 (Right). Notice in the figure very small undulations in amplitude. 

To construct snake solitons with amplitudes large enough, we had to move away from the Hopf curve, 
as much as possible, but at the same time to be sure not to be outside of the parameters ranges for the 
existence of the snaking soliton. That could be achieved by varying one or more of the system parameters. 
First, we varied the first bifurcation parameter e, slowly away from the Hopf curve where 63 = b^Hopf, and 
e = tHopf, and we noticed that the snaking soliton had very small amplitude of Ai{t). Since the soliton in 
this case had only a magnitude of only 10^^, we decided to vary the second bifurcation parameter 63, which 
stands for the cubic gain when negative. We found that the domain of existence for the snaking soliton as a 
function of 63 was [—2.25234,-0.143456], passing through the Hopf curve value of &3Hop/ = —1.89646. For 
the unperturbed values of 63 — —0.835, and e = —0.1, we numerically integrate the Euler-Lagrange system 
of ODEs (11.17) using the Matlab's built-in function ode45 and we plot the periodic orbit, which is shown in 
Fig. 6 (Left), while the resulting periodic time series for Ai{t), a{ t), and (/)(i) are shown in Fig. 6 (Right). 
The resulting time series are then substituted back in the ansatz ( |II.13 ) and the spatiotemporal structure 
\A{x^ t) \ of the snaking soliton obtained by the variational formulation is presented in Fig. 7. As the various 
system parameters ci, C3, C5, hi, 63, 65, and e within the stable regime are varied, the effects of the amplitude, 
position, width (and, less importantly, phase) may also be studied, and this is discussed subsequently. As 
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FIG. 6: Numerical simulation using the variational approximation for the snake of Fig. 1: periodic orbit or limit 
cycle (Left), periodic time series for Ai{t), a{t), and 0(t) (Right) 




FIG. 7: Snake soliton using the ODEs (variational approximation): unperturbed case 



observed from the simulations of the full PDE and ODEs system, the snaking soliton obtained from the 
ODE system (11.17 1 qualitatively agree for a short time to the solution of the complex nonlinear PDE (II. 1 1. 
Let us take a closer look to the snake solitons shown in Figs. 7-10. The soliton now "snakes" or wiggles as 



its position varies periodically in time. Note that the amplitude \A{t)\ varies periodically as Ai{t) varies'; 
but there would be additional amplitude modulation due to the periodic variation of (j){t). Next, we shall 
consider the effect of all the various parameters in the CCQGLE (II. 1) on the shape (amplitude, position, 
phase, period) and stability of the snake. 

In considering the parameter effects on snake shape and period, note that the wave is a spatially coherent 
structure (or a "collective coordinate" given by the trial function) whose parameters oscillate in time. Hence, 
the temporal period of the snake is the same as the period T of the oscillations of Ai{t), (f>{t), and a{t) on 
their limit cycle. As for the peak amplitude and peak position of the snake, these are determined by the 
peak amplitude Aip of Ai{t), and the peak position (pp of (/)(t) respectively. Notice that from (11.13) we 
can regard the width and the amplitude of the snake as being inversely proportional to position (p{t) for the 
snake i.e., at any time t when the amplitude is minimum, the width will be minimum, so the position is 
maximum and vice versa. So, maximum deflection from the horizontal position x = const, is obtained when 
the position of the snake is maximum, and hence the width and amplitude are minimum. 

Keeping the above in mind, we vary the parameters of the CCQGLE in turn and we observe the resulting 
effects on Aip (the peak amplitude), (j)p (the position), and T (the temporal period) of the snake soliton: 



(i) For increased 61, the values of A 



and T all increase. 



(ii) Increasing 63 augments all of ^ip, (j^p, and T. 



(iii) Increasing 65 increases all of Aip, (f>p, and T. 



(iv) Raising C\ increases Aip, (pp, but decreases T. 



(v) Incrementing C3 decreases all of Aip, (f>p, and T. 



(vi) Augmenting C5 causes a decrease in Aip, cpp, and increases T. 



(vii) Raising e causes Aip, Op to rise, but T to fall. 





FIG. 8: Snake soliton using the ODEs (variational approximation): Left 63 = —0.6, Right 65 = 0.62 



The above constitute our detailed predictions of the various parameters in the CCQGLE on the amplitude, 
position, and temporal width of the snake solitons. We have verified that each set of predictions (i)-(vii) 
above agree when the corresponding parameter is varied in the solitary wave simulation for the full PDE. 
Note also that Aiit) and (j){t) are always in phase, so that Aip and (pp occur simultaneously. Thus, the 



snaking solitons are tallest where they have the greatest width, 
simulation, as well as those in [3]. 



This is completely consistent with our 
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FIG. 9: Snake soliton using the ODEs (variational approximation): Left ci = 0.55, Right C3 = 2 




FIG. 10: Snake sohton using the ODEs (variational approximation): Left C5 = 0.8, Right e — —0.08 



VI. CONCLUSION 

In this paper we discussed a brief theoretical framework for analyzing the full spatiotemporal structure of 
one class of solitary waves (snaking soliton) in the CCQGLE and presented some numerical simulations of 
dissipative snaking soliton solutions to the CCQGLE. The results obtained analytically using the variational 
approximation for the snaking soliton is compared with the numerical simulations of the CCQGLE. The 
specific theoretical modeling includes the use of a recent variational formulation and significantly generalized 
trial function for the solitary waves solutions. In addition, the resulting Euler-Lagrange equations are 
treated in an entirely different way by looking at their stable periodic solutions (or limit cycles) resulting 
from supercritical Hopf bifurcations. Oscillations of their trial function parameters on these limit cycles 
provide the pulsations of the amplitude, width/position, and phase of the solitons. The model also allows 
for detailed predictions regarding the other types of solitons i.e., pulsating, chaotical and creeping. 

In the future, we will explore on the (2D) solutions called spinning solitons as well as the (3D) solutions 
called optical bullets. These depict a confined spatiotemporal soliton in which the balance between the 
focusing nonlinearity and the spreading while propagating through medium provides the shape of a bullet 

mm- 
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